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(4)  Introduction 


Previous  studies  have  found  that  there  is  a  strong  correlation  between  mamnaographic 
breast  density  and  the  risk  of  breast  cancer.  Mammographic  breast  density  has  been  used  by 
researchers  in  many  studies  to  estimate  breast  cancer  risk  of  epidemiological  factors,  monitor  the 
effects  of  preventive  treatments  such  as  tamoxifen  or  dietary  interventions,  monitor  the  breast 
cancer  risk  of  hormone  replacement  therapy,  and  investigate  factors  affecting  mammographic 
sensitivity  and  cancer  prognosis.  However,  most  studies  used  Breast  Imaging  Reporting  and 
Data  System  (BI-RADS)  density  rating  as  a  measure  of  mammographic  breast  density,  which 
contributes  large  inter-  and  intraobserver  variations  and  may  reduce  the  sensitivity  of  the 
analysis. 


The  goal  of  this  proposed  project  is  to  develop  a  fully  automated  technique  to  assist 
radiologists  in  estimating  mammographic  breast  density.  We  hypothesize  that  the  computerized 
technique  can  accurately  and  efficiently  segment  the  dense  area  on  digitized  or  digital 
mammograms,  thereby  eliminating  inter-  and  intra-observer  variations.  The  dense  area  as  a 
percentage  of  total  breast  area  thus  estimated  will  be  more  consistent  and  reproducible  than 
radiologists'  subjective  BI-RADS  rating.  To  accomplish  this  goal,  we  will  (1)  collect  a  large 
database  of  normal  mammograms,  including  digitized  film  mammograms  and  digital 
mammograms,  for  training  and  testing  the  dense  area  segmentation  program;  (2)  evaluate  the 
correlation  between  the  radiologists’  breast  density  classification  based  on  BI-RADS  lexicon  and 
the  percent  breast  dense  area;  and  (3)  study  the  correlation  of  percent  breast  dense  area  between 
different  views  of  the  same  breast  and  between  the  same  view  of  the  two  breasts;  and  (4) 
investigate  the  correlation  between  the  percent  breast  dense  area  estimated  from  mammograms 
and  the  volumetric  dense  breast  tissue  estimated  from  a  data  set  of  magnetic  resonance  (MR) 
breast  images.  These  comparisons  will  provide  important  information  regarding  the  consistency 
of  the  BI-RADS  rating  with  the  measured  percent  breast  dense  area,  the  appropriate  measure  of 
%  dense  area  from  different  mammographic  views,  and  the  usefulness  of  using  the  percent  breast 
dense  area  on  mammograms  as  an  indicator  of  volumetric  breast  tissue  density. 

It  is  expected  that  this  project  will  produce  a  fully  automated  and  effective  tool  for  analysis 
of  mammographic  breast  density,  which  can  be  applied  to  routinely  acquired  mammograms 
without  special  calibrations.  This  will  facilitate  studies  of  various  factors  associated  with  breast 
cancer  risk  and  mammographic  sensitivity,  and  monitoring  the  effects  of  interventional  or 
preventive  strategies.  The  image  analysis  tool  will  therefore  contribute  to  the  understanding  of  the 
relationship  of  density  to  breast  cancer  risk,  detection,  prognosis,  and  to  the  prevention  and 
treatment  of  breast  cancer. 
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(5)  Body 

In  the  current  project  year  (7/1/01-6/30/02),  we  have  performed  the  following  studies: 

(A)  Collection  of  Database 

We  have  obtained  the  research  log  (with  IRB  approval)  of  the  previous  research  project 
that  developed  a  dynamic  MR  imaging  technique  for  the  characterization  of  breast  lesions. 
Based  on  the  record,  we  found  the  mammographic  files  of  69  patients  who  have  corresponding 
MR  breast  images  and  mammograms  that  are  useful  for  the  current  study.  The  other  cases  have 
incomplete  records,  such  as  mammograms  missing  from  the  patient  files,  certain  MR  sequences 
not  available.  We  will  continue  to  track  these  cases  to  see  if  some  of  the  missing  information  can 
be  located.  We  have  developed  a  database  to  categorize  the  cases  for  our  project.  The  database 
includes  the  digitized  mammograms,  the  corresponding  MR  images,  and  the  pathological 
information  obtained  from  biopsy.  We  have  digitized  the  mammograms  and  archived  them  in 
our  jukebox.  For  the  MR  study,  we  have  retrieved  the  digital  files  from  the  archive  in  the  MR 
laboratory,  transferred  all  the  files  to  our  jukebox,  and  extracted  the  pre-contrast  sequences  for 
our  analysis. 

The  MR  image  series  used  in  this  study  included  coronal  3D  Ti-weighted  pre-contrast 
images  acquired  with  the  following  technique:  coronal  sections  2-5  mm  thick,  32  slices;  3D 
Spoiled  Gradient-Recalled  Echo  (SPGR);  TE=3.3ms;  TR=10ms,  Flip=40°,  matrix=256xl28, 
FOV=28-32cm  right/left,  14-16cm  superior/inferior,  scan  time=2  min  38  sec.  This  3D  SPGR 
sequence  produced  full  volume  coverage  of  both  breasts  with  contiguous  image  slices.  An 
example  of  images  from  one  breast  is  shown  in  Fig.  1.  Although  this  is  not  the  optimal  pulse 
sequence  for  separating  water  and  fat,  the  fibroglandular  parenchyma  (-water)  and  fatty  tissue  are 
well  separated  with  the  heavily  Ti-weighted  acquisition  and  therefore  the  series  was  chosen  for 
the  study. 

(B)  Development  of  Graphical  User  Interface  (GUI)  for  Analysis  of  MR  Images 

We  have  developed  a  GUI  for  analysis  of  the  MR  images  semi-automatically.  The  GUI 
was  written  based  on  Matlab  in  an  AlphaStation.  The  GUI  first  displays  all  slices  of  the  pre¬ 
contrast  MR  sequence  of  one  breast  and  the  corresponding  mammogram  of  the  same  breast  (Fig. 
1).  The  user  can  switch  between  the  craniocaudal  (CC)  or  the  mediolateral  oblique  (MLO)  view 
mammogram  for  comparison.  The  computer  then  displays  each  MR  slice  sequentially  (Fig.  2). 
For  each  slice  breast  boundary  is  first  detected  automatically.  A  low  pass  filtering  is  employed 
within  the  breast  boundary  to  estimate  the  low-frequency  gray  level  background  trend  due  to  the 
field  non-uniformity  of  the  breast  coil.  The  low  frequency  background  is  subtracted  from  the 
original  image.  After  background  correction,  the  histogram  within  the  breast  region  is 
determined  and  interactive  thresholding  is  used  to  separate  the  fibroglandular  tissue  region 
(displayed  as  dark)  from  the  fatty  tissue  (white)  region.  After  the  segmentation  is  performed  for 
all  slices,  the  total  fibroglandular  tissue  volume  in  the  breast  is  estimated  by  summing  the  dark 
voxels  over  all  slices,  and  the  total  breast  volume  is  estimated  by  summing  all  the  voxels  within 
the  breast.  The  percent  fibroglandular  tissue  volume  can  then  be  derived  as  the  ratio  of  the  total 
fibroglandular  tissue  volume  to  the  total  breast  volume.  Some  of  the  segmentation  steps  are 
detailed  below: 

Breast  segmentation 
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We  performed  the  analysis  for  each  individual  breast.  The  breast  region  in  the  MR  slices 
of  the  anterior  portion  of  the  breast  can  be  easily  segmented  from  the  background  by  automatic 
edge  detection.  As  the  slices  get  close  to  the  chest  wall,  the  breast  tissue  merges  with  the  chest 
tissue.  Often  there  is  no  distinguishable  boundary  between  the  breast  tissue  and  the  chest  tissue 
in  these  slices.  An  expert  breast  radiologist  will  inspect  the  computer  segmentation  on  the  GUI 
slice  by  slice  and  correct  the  edge  location  if  necessary  (Fig.  3).  In  our  GUI,  several  options  for 
boundary  correction  are  available.  The  first  option  is  to  use  a  snake  algorithm  to  push  the  initial 
boundary  to  be  closer  to  the  actual  boundary.  The  user  can  determine  when  to  stop  the  iteration. 
The  second  option  is  to  use  an  active  contour  algorithm  to  push  the  initial  boundary  to  the  actual 
boundary.  This  is  similar  to  the  snake  algorithm  except  that  the  energy  terms  are  different.  The 
third  option  is  to  manually  move  the  edge,  which  is  represented  by  a  polygon  with  an  edge  point 
at  each  apex,  to  the  user  chosen  location.  The  last  option  is  to  manually  draw  a  new  edge  by 
choosing  the  apex  points  of  a  polygon. 

Background  correction 

Background  correction  is  applied  to  the  breast  region  within  the  detected  boundary.  We 
have  evaluated  two  background  correction  methods:  (1)  A  low  pass  filtering  method  in  which 
the  low  frequency  gray  level  background  of  the  breast  is  estimated  by  low-pass  filtering  and 
subtracted  from  the  original  image.  (2)  A  distance- weighted  background  estimation  method  in 
which  a  given  pixel  in  the  background  image  is  calculated  as  the  weighted  sum  of  the  four 
boundary  pixels.  The  weights  are  estimated  by  the  inverse  distance  from  the  breast  boundary  to 
the  given  pixel.  The  background  image  is  then  subtracted  from  the  original  image. 

After  comparison  with  a  large  number  of  images,  it  was  found  that  the  second  method  is 
more  robust  than  the  first.  The  first  method  tends  to  overestimate  the  background  for  images 
containing  a  large  area  of  fibroglandular  tissue.  The  second  method  can  closely  approximate  the 
background  gray  level  trend  and  correct  for  the  non-uniform  background  in  the  original  image. 

Gray  level  thresholding 

Gray  level  thresholding  is  used  to  separate  the  fibroglandular  tissue  from  the  fatty  tissue. 
We  employ  an  interactive  thresholding  method  in  which  the  gray  level  histogram  of  the  breast 
region  is  displayed  on  the  GUI.  The  thresholded  image  where  the  dense  tissue  is  displayed  as 
dark  voxels  and  fatty  tissue  is  displayed  in  white  is  also  displayed  and  updated  instantly 
whenever  the  threshold  is  changed.  By  comparing  the  thresholded  image  and  the  original  image 
side-by-side,  the  user  can  adjust  the  threshold,  using  the  threshold  slider  under  the  histogram, 
until  the  fibroglandular  tissue  is  extracted  from  the  fatty  background. 

Density  mask 

In  many  MR  slices,  the  skin  forms  a  ring  of  dark  voxels  around  the  breast  boundary.  To 
exclude  these  skin  voxels  from  being  counted  as  fibroglandular  tissue,  a  density  mask  is  defined 
by  applying  morphological  erosion  to  the  breast  boundary  to  obtain  a  mask  slightly  smaller  than 
the  breast  region.  The  user  can  choose  the  structure  element  size  for  the  morphological 
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operation.  The  resulting  density  mask  is  displayed  superimposed  on  the  thresholded  breast  image 
so  that  the  user  can  judge  if  the  skin  voxels  are  excluded  from  the  mask  region. 

Percent  fibroglandular  tissue  volume  calculation 


After  the  segmentation  process  is  applied  to  all  slices  of  the  MR  sequence  of  the  breast, 
all  voxels  within  the  density  mask  in  each  slice  are  summed  over  the  entire  breast  to  obtain  the 
fibroglandular  tissue  volume  in  the  breast.  Similarly,  the  voxels  within  the  breast  boundary  in 
each  slice  are  summed  over  all  slices  to  obtain  the  total  volume  of  the  breast.  The  ratio  of  the 
fibroglandular  tissue  volume  to  the  total  breast  volume  yields  the  percent  volumetric 
fibroglandular  tissue  in  the  breast. 


File  Edit  Tools  Window  Help 
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Fig.  1.  The  first  screen  of  the  graphical  user  interface  that  we  developed  for  analysis  of  the  MR 
breast  images.  The  analysis  was  performed  one  breast  at  a  time.  The  displayed  slices  on 
the  left  showed  the  MR  sequence  from  the  nipple  to  the  chest  wall  of  a  coronal  MR 
image  sequence.  The  corresponding  mammogram  of  the  same  breast  is  shown  in  the 
lower  middle.  The  mammogram  can  be  switched  between  the  CC  view  and  the  MLO 
view. 
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Fig.  2.  The  second  screen  of  the  GUI.  The  MR  sequence  is  analyzed  one  slice  at  a  time.  The 
breast  boundary  was  detected.  The  non-uniform  gray  level  background  was  estimated  by 
low-pass  filtering  and  subtracted  from  the  original  image.  The  background  corrected 
image  is  shown  in  the  middle.  The  histogram  of  this  image  is  obtained  and  displayed  at 
the  bottom  of  the  screen.  The  user  interactively  thresholds  the  image.  The  image  results 
from  thresholding  is  shown  on  the  right.  The  user  subjectively  determines  the  proper 
threshold  for  separately  the  fibroglandular  tissue  from  the  fatty  tissue. 
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Fig.  3.  The  third  screen  of  the  GUI.  If  the  user  is 
not  satisfied  with  the  breast  boundary 
displayed  on  the  second  screen,  the  user 
can  click  a  boundary  correction  button  on 
the  second  screen  and  get  to  this  third 
screen.  In  this  screen,  the  automatically 
detected  edge  is  displayed  on  the  breast 
slice.  The  user  can  use  a  snake  algorithm 
to  move  the  edge  automatically,  or 
manually  adjust  the  edge  location  by 
moving  the  edge  points  individually. 


Fig.  4.  The  fourth  screen  of  the  GUI.  After 
edge  detection,  background  correction, 
and  thresholding,  the  user  can  generate  a 
“density  mask”  to  define  the  region 
where  the  thresholded  black  voxels 
representing  the  fibroglandular  tissue 
should  be  summed.  The  mask  excludes 
the  skin  layer  from  the  black  voxels  by 
applying  morphological  erosion  to  the 
breast  boundary. 


(C)  Estimation  of  Percent  Dense  Area  on  Mammograms 

We  are  developing  an  automated  image  analysis  tool.,  Mammographic  Density 
ESTimator  (MDEST),  to  assist  radiologists  in  estimating  mammographic  breast  density. 
MDEST  performs  dynamic  range  compression,  breast  boundary  tracking,  pectoral  muscle 
segmentation  for  the  MLO  view,  automatic  thresholding  based  on  gray  level  histogram  analysis, 
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and  calculates  the  percent  dense  area  on  a  mammogram.  In  a  preliminary  study,  we  found  that 
the  correlation  between  the  computer-estimated  percent  dense  area  and  radiologists’  manual 
segmentation  was  0.94  and  0.91,  respectively,  for  CC  and  MLO  views,  with  a  mean  bias  of  less 
than  2%.  As  a  result  of  this  study,  a  journal  article  was  published  in  the  Medical  Physics  Journal 
[ref.  1]. 


(D)  Correlation  between  Percent  Dense  Area  and  Percent  Volumetric  Fibroglandular 
Tissue 

To  date,  we  have  performed  segmentation  of  the  MR  images  and  their  corresponding 
mammograms  for  37  cases  [ref.  2].  Scatter  plots  of  the  percent  volumetric  fibroglandular  tissue 
versus  the  percent  dense  area  on  mammograms  are  shown  in  Fig.  5(a)  and  5(b)  for  the  CC-  and 
MLO-view  mammograms,  respectively.  The  correlation  of  percent  dense  area  of  the  CC  and 
MLO  views  with  the  percent  volumetric  fibroglandular  tissue  on  MR  images  was  found  to  be 
0.93  and  0.91,  respectively,  with  a  mean  bias  of  4.4%. 


Fig.  5.  Correlation  of  %  volumetric  fibroglandular  tissue  on  MR  images  with  %  dense  area  on 
mammograms  for  37  patients.  The  left  and  right  breasts  are  plotted  as  separate  data  points 
on  each  graph.  The  dash  lines  are  linear  least  squares  fits  to  the  data  points. 


(6)  Key  Research  Accomplishments 

•  Collect  the  database  of  MR  breast  image  sequence  and  the  corresponding  mammograms  from 
the  previous  MR  dynamic  imaging  study  (Task  1  and  Task  3) 

•  Develop  the  Graphical  User  Interface  (GUI)  for  interactive  segmentation  of  the  MR  images 
(Task  3). 

•  Develop  preliminary  methods  for  automated  estimation  of  mammographic  density  (Task  2). 
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•  Analyze  37  cases  of  the  MR  images  and  mammograms  and  study  the  correlation  between  the 
percent  volumetric  fibroglandular  tissue  and  percent  mammographic  dense  area  (Task  3). 

•  Present  preliminary  results  at  the  International  Workshop  for  Digital  Mammography  (IWDM) 
in  Bremen,  Germany. 


(7)  Reportable  Outcomes 

As  a  result  of  the  support  by  the  USAMRMC  BCRP  grant,  we  have  conducted  studies  to 
investigate  the  correlation  between  mammographic  density  and  MR  volumetric  fibroglandular 
tissue.  The  preliminary  results  were  presented  in  an  international  conference.  A  conference 
proceeding  paper  will  be  published  by  the  IWDM.  In  addition,  we  have  developed  a  preliminary 
method  for  automated  estimation  of  mammographic  density,  and  performed  a  study  to  evaluate 
the  correlation  between  the  computer-determined  mammographic  density  and  radiologists’ 
manual  segmentation  using  interactive  thresholding.  The  results  were  published  in  an  article  in 
the  Medical  Physics  Journal. 

Journal  Articles: 


Zhou  C,  Chan  HP,  Petrick  N,  Helvie  MA,  Goodsitt  MM,  Sahiner  B,  Hadjiiski  LM. 
Computerized  image  analysis:  Estimation  of  breast  density  on  mammograms.  Medical  Physics 
2001;28:1056-1069. 

Conference  Proceeding: 

Chan  HP,  Hadjiiski  LM,  Roubidoux  MA,  Helvie  MA,  Paquerault  S,  Sahiner  B,  Chenevert  T, 
Goodsitt  MM.  Breast  density  estimation:  correlation  of  mammographic  density  and  MR 
volumetric  density.  Proceedings  of  the  6th  International  Workshop  on  Digital  Mammography. 
IWDM-2002.  Bremen,  Germany.  June  22-25,  2002.  (in  press). 

Conference  Presentation: 


Chan  HP,  Hadjiiski  LM,  Roubidoux  MA,  Helvie  MA,  Paquerault  S,  Sahiner  B,  Chenevert  T, 
Goodsitt  MM.  Breast  density  estimation:  correlation  of  mammographic  density  and  MR 
volumetric  density.  Poster  presentation  at  the  6th  International  Workshop  on  Digital 
Mammography.  IWDM-2002.  Bremen,  Germany.  June  22-25, 2002. 


(8)  Conclusions 

During  this  project  year,  we  have  collected  the  database  of  cases  that  have  corresponding 
magnetic  resonance  (MR)  images  and  mammograms  acquired  in  a  previous  project.  We  have 
developed  methods  and  a  graphical  user  interface  to  analyze.  MR  breast  images.  The 
fibroglandular  tissue  in  the  MR  breast  images  was  segmented  with  a  semi-automatic  method. 
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We  are  also  developing  methods  and  computer  programs,  referred  to  as  Mammographic  Density 
ESTimator  (MDEST),  to  automatically  segment  digitized  mammograms  and  estimate  the  percent 
dense  area  on  CC  and  MLO  views.  The  preliminary  results  for  analyzing  37  cases  indicate 
strong  correlation  between  mammographic  density  and  the  volumetric  fibroglandular  tissue 
estimated  from  MR  breast  images. 

These  results,  although  still  preliminary,  indicate  that  the  usefulness  of  mammographic 
density  as  a  surrogate  for  breast  density  estimation.  The  automated  analysis  will  provide  a 
reproducible  and  efficient  method  for  analysis  of  mammographic  density.  When  fully  developed, 
MDEST  is  expected  to  contribute  to  the  understanding  of  the  relationship  of  mammographic 
density  to  breast  cancer  risk,  detection,  and  prognosis,  and  to  the  prevention  and  treatment  of 
breast  cancer. 
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Computerized  image  analysis:  Estimation  of  breast  density 
on  mammograms 

Chuan  Zhou,  Heang-Ping  Chan,a)  Nicholas  Petrick,  Mark  A.  Helvie,  Mitchell  M.  Goodsitt, 
Berkman  Sahiner,  and  Lubomir  M.  Hadjiiski 

Department  of  Radiology,  The  University  of  Michigan,  Ann  Arbor,  Michigan  48109-0030 
(Received  15  September  2000;  accepted  for  publication  4  April  2001) 

An  automated  image  analysis  tool  is  being  developed  for  the  estimation  of  mammographic  breast 
density.  This  tool  may  be  useful  for  risk  estimation  or  for  monitoring  breast  density  change  in 
prevention  or  intervention  programs.  In  this  preliminary  study,  a  data  set  of  4- view  mammograms 
from  65  patients  was  used  to  evaluate  our  approach.  Breast  density  analysis  was  performed  on  the 
digitized  mammograms  in  three  stages.  First,  the  breast  region  was  segmented  from  the  surrounding 
background  by  an  automated  breast  boundary-tracking  algorithm.  Second,  an  adaptive  dynamic 
range  compression  technique  was  applied  to  the  breast  image  to  reduce  the  range  of  the  gray  level 
distribution  in  the  low  frequency  background  and  to  enhance  the  differences  in  the  characteristic 
features  of  the  gray  level  histogram  for  breasts  of  different  densities.  Third,  rule-based  classification 
was  used  to  classify  the  breast  images  into  four  classes  according  to  the  characteristic  features  of 
their  gray  level  histogram.  For  each  image,  a  gray  level  threshold  was  automatically  determined  to 
segment  the  dense  tissue  from  the  breast  region.  The  area  of  segmented  dense  tissue  as  a  percentage 
of  the  breast  area  was  then  estimated.  To  evaluate  the  performance  of  the  algorithm,  the  computer 
segmentation  results  were  compared  to  manual  segmentation  with  interactive  thresholding  by  five 
radiologists.  A  “true”  percent  dense  area  for  each  mammogram  was  obtained  by  averaging  the 
manually  segmented  areas  of  the  radiologists.  We  found  that  the  histograms  of  6%  (8  CC  and  8 
MLO  views)  of  the  breast  regions  were  misclassified  by  the  computer,  resulting  in  poor  segmen¬ 
tation  of  the  dense  region.  For  the  images  with  correct  classification,  the  correlation  between  the 
computer-estimated  percent  dense  area  and  the  “truth”  was  0.94  and  0.91,  respectively,  for  CC  and 
MLO  views,  with  a  mean  bias  of  less  than  2%.  The  mean  biases  of  the  five  radiologists’  visual 
estimates  for  the  same  images  ranged  from  0.1%  to  11%.  The  results  demonstrate  the  feasibility  of 
estimating  mammographic  breast  density  using  computer  vision  techniques  and  its  potential  to 
improve  the  accuracy  and  reproducibility  of  breast  density  estimation  in  comparison  with  the 
subjective  visual  assessment  by  radiologists.  ©  2001  American  Association  of  Physicists  in  Medi¬ 
cine .  [DOI:  10.1118/1.1376640] 

Key  words:  mammography,  computer-aided  diagnosis,  breast  density,  breast  cancer  risk,  image 
segmentation,  thresholding 


I.  INTRODUCTION 

Breast  cancer  is  one  of  the  leading  causes  for  cancer  mortal¬ 
ity  among  women.1  One  in  every  eight  women  will  develop 
breast  cancer  at  some  point  in  their  lives.  The  most  success¬ 
ful  method  for  the  early  detection  of  breast  cancer  is  screen¬ 
ing  mammography.  Currently,  mammograms  are  analyzed 
visually  by  radiologists.  Because  of  the  subjective  nature  of 
visual  analysis,  qualitative  responses  may  vary  from  radiolo¬ 
gist  to  radiologist.  Therefore,  a  computerized  method  for 
analyzing  mammographic  features  would  be  useful  as  a 
supplement  to  the  radiologist’s  assessment.  Previous  re¬ 
search  efforts  in  computer-aided  diagnosis  (CAD)  for  breast 
cancer  detection  mainly  concentrated  on  detection  and  char¬ 
acterization  of  masses  and  microcalcifications  on  mammo¬ 
grams  by  using  computer  vision  techniques.  It  has  been  dem¬ 
onstrated  that  an  effective  CAD  algorithm  can  improve  the 
diagnostic  accuracy  of  breast  cancer  characterization  on 
mammograms,  which,  in  turn,  may  reduce  unnecessary  biop¬ 
sies.  In  this  work,  we  are  studying  the  feasibility  of  develop¬ 


ing  a  CAD  system  for  an  analysis  of  breast  density  on  mam¬ 
mograms.  Studies  have  shown  that  there  is  a  strong  positive 
correlation  between  breast  parenchymal  density  on  mammo¬ 
grams  and  breast  cancer  risk.2-9  The  relative  risk  is  estimated 
to  be  about  4  to  6  times  higher  for  women  whose  mammo¬ 
grams  have  parenchymal  densities  over  60%  of  the  breast 
area,  as  compared  to  women  with  less  than  5%  of  parenchy¬ 
mal  densities. 

An  important  difference  between  breast  density  as  a  risk 
factor  and  most  other  risk  factors  is  the  fact  that  breast  tissue 
density  can  be  changed  by  dietary  or  hormonal 
interventions.6,10,11  Although  there  is  no  direct  evidence  that 
changes  in  mammographic  breast  densities  will  lead  to 
changes  in  breast  cancer  risk,  the  strong  correlation  between 
breast  density  and  breast  cancer  risk  has  prompted  research¬ 
ers  to  use  mammographic  density  as  an  indicator  for  moni¬ 
toring  the  effects  of  intervention  as  well  as  for  studying 
breast  cancer  etiology.6,11-13 

Different  methods  have  been  used  for  the  evaluation  of 
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mammographic  breast  density.  Earlier  studies  used  a  subjec¬ 
tive  visual  assessment  of  the  breast  parenchyma  primarily 
based  on  the  four  patterns  described  by  Wolfe2  (N1  is  com¬ 
prised  entirely  of  fat;  PI  has  up  to  25%  nodular  densities;  P2 
has  over  25%  nodular  mammographic  densities;  DY  contains 
extensive  regions  of  homogeneous  mammographic  densi¬ 
ties).  The  subjectivity  in  classifying  the  mammographic  pat¬ 
terns  introduced  large  variability  in  the  risk  estimation.  Later 
studies  used  more  quantitative  estimates,  such  as  planimetry, 
to  measure  the  dense  area  in  the  breast  manually  outlined  by 
radiologists  on  mammograms.3’7  These  studies  indicate  that 
the  percentage  (%)  of  mammographic  densities  relative  to 
the  breast  area  can  predict  the  breast  cancer  risk  more  accu¬ 
rately  than  a  qualitative  assessment  of  mammographic  pat¬ 
terns.  Warner  et  al 15  conducted  a  meta-analysis  of  the  stud¬ 
ies  published  between  1976  and  1990  to  investigate  the 
effect  of  different  methods  of  classification  on  estimates  of 
cancer  risk.  They  found  that  the  mammographic  parenchy¬ 
mal  pattern  does  correlate  with  the  breast  cancer  risk.  The 
magnitude  of  the  risk  varies  according  to  the  method  used  to 
evaluate  the  mammograms.  With  the  quantitative  estimates 
of  mammographic  density,  the  difference  in  risk  between  the 
highest  and  the  lowest  risk  category  is  substantial  and  is 
greater  than  the  risks  associated  with  most  other  risk  factors 
for  breast  cancer.  More  recent  studies  used  fractal  texture 
and  the  shape  of  the  gray  level  histogram14  to  quantify  the 
parenchymal  pattern  or  used  interactive  thresholding  on  digi¬ 
tized  mammograms  to  segment  the  dense  area.11,15  It  was 
reported  that  the  thresholding  method  provided  a  higher  risk 
value  than  the  texture  measure  or  the  histogram  shape.16 
Other  researchers  have  attempted  to  calculate  a  breast  den¬ 
sity  index  to  model  the  radiologists’  perception.17 

In  clinical  practice,  radiologists  routinely  estimate  the 
breast  density  on  mammograms  by  using  the  BI-RADS  lexi¬ 
con  as  recommended  by  the  American  College  of 
Radiology18  in  order  to  provide  a  reference  for  mammo¬ 
graphic  sensitivity.  Because  of  the  lack  of  a  quantitative 
method  for  breast  density  estimation,  researchers  often  use 
the  BI-RADS  rating  for  monitoring  responses  to  preventive 
or  interventional  treatment  and  the  associated  changes  in 
breast  cancer  risk.19  We  have  found  that  there  is  a  large 
interobserver  variability  in  the  BI-RADS  ratings  among  ex¬ 
perienced  mammographers. 20,21  An  automated  and  quantita¬ 
tive  estimation,  as  investigated  in  this  study,  will  provide  not 
only  an  efficient  means  to  measure  mammographic  density, 
but  also  a  reproducible  estimate  that  will  reduce  the  inter- 
and  intraobserver  variability  of  mammographic  density  mea¬ 
surements.  This  image  analysis  tool  will  therefore  allow  re¬ 
searchers  to  study  more  definitively  the  relationship  of  mam¬ 
mographic  density  to  breast  cancer  risk,  detection,  prognosis, 
and  mammographic  sensitivity,  and  to  better  monitor  the  re¬ 
sponse  of  a  patient  to  preventive  or  interventional  treatment 
of  breast  cancers. 

In  this  paper,  we  will  describe  the  image  processing  tech¬ 
niques  used  in  our  automated  breast  density  segmentation 
algorithm.  The  performance  of  the  computer  segmentation 
was  evaluated  by  a  comparison  with  the  average  segmenta- 
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tion  by  5  radiologists  using  interactive  thresholding  in  the 
same  data  set. 

II.  MATERIALS  AND  METHODS 

A.  Database 

A  data  set  consisting  of  260  mammograms  of  65  patients 
was  used  for  the  development  of  the  histogram  analysis 
method  in  this  study.  Each  case  contains  the  craniocaudal 
(CC)  view  and  the  mediolateral  oblique  (MLO)  view  of  both 
breasts  of  the  patient.  The  first  50  mammograms  were  con¬ 
secutive  screening  cases  from  the  patient  files  in  the  Radiol¬ 
ogy  Department  at  the  University  of  Michigan.  After  data 
analysis,  it  was  found  that  there  were  very  few  dense  breasts 
in  the  initial  data  set.  An  additional  15  cases  visually  judged 
by  radiologists  to  be  dense  breasts  were  then  randomly  se¬ 
lected  and  mixed  with  the  initial  set.  The  images  were  pro¬ 
cessed  individually  without  knowing  their  BI-RADS  catego¬ 
ries.  The  mammograms  were  acquired  with  mammography 
systems  approved  by  the  Mammography  Quality  Standards 
Act  (MQSA)  and  were  digitized  with  a  LUMISYS  85  laser 
film  scanner  with  a  pixel  size  of  50/rmX50/im  and  4096 
gray  levels.  The  gray  levels  are  linearly  proportional  to  op¬ 
tical  densities  (O.D.)  from  0.1  to  greater  than  3  O.D.  units. 
The  nominal  O.D.  range  of  the  scanner  is  0-4  with  large 
pixel  values  in  the  digitized  mammograms  corresponding  to 
low  O.D.  The  full  resolution  mammograms  were  first 
smoothed  with  a  16X 16  box  filter  and  subsampled  by  a  fac¬ 
tor  of  16,  resulting  in  800//,mX800/zm  images  of  approxi¬ 
mately  225X300  pixels  in  size  for  small  films  and  300 
X375  pixels  for  large  films. 

B.  Breast  segmentation  and  image  enhancement 

The  breast  image  is  first  segmented  from  the  surrounding 
image  background  by  boundary  detection.  The  detected 
boundary  separated  the  breast  from  other  background  fea¬ 
tures  such  as  the  directly  exposed  area,  patient  identification 
information,  and  lead  markers.  The  density  analysis  was  per¬ 
formed  only  within  the  breast  region.  An  automated  breast 
boundary  tracking  technique  developed  previously22,23  was 
modified  to  improve  its  performance.  Briefly,  the  technique 
used  a  gradient-based  method  to  search  for  the  breast  bound¬ 
ary.  The  background  of  the  image  was  estimated  initially  by 
searching  for  the  largest  background  peak  from  the  gray 
level  histogram  of  the  image.  After  subtracting  this  back¬ 
ground  level  from  the  breast  region,  a  simple  edge  was  found 
by  a  line-by-line  gradient  analysis  from  the  top  to  the  bottom 
of  the  image.  The  criterion  used  in  detecting  the  edge  points 
was  the  steepness  of  the  gradient  of  four  adjacent  pixels 
along  the  horizontal  direction.  The  steeper  the  gradient,  the 
greater  the  likelihood  that  an  edge  existed  at  that  correspond¬ 
ing  image  point.  The  simple  edge  served  as  a  starting  point 
for  a  more  accurate  tracking  algorithm  that  followed.  The 
tracking  of  the  breast  boundary  started  from  approximately 
the  middle  of  the  breast  image  and  moved  upward  and  down¬ 
ward  along  the  boundary.  The  direction  to  search  for  a  new 
edge  point  was  guided  by  the  previous  edge  points.  The  edge 


1058 


Zhou  et  a/.:  Computerized  image  analysis 


1058 


Fig.  1.  (a)  A  mammogram  from  our 
image  database;  (b)  the  image  super¬ 
imposed  with  the  detected  breast 
boundary  and  pectoral  muscle  bound¬ 
ary;  (c)  the  binary  map  of  the  seg¬ 
mented  breast  region. 


location  was  again  determined  by  searching  for  the  maxi¬ 
mum  gradient  along  the  gray  level  profile  normal  to  the 
tracking  direction.  Since  the  boundary  tracking  was  guided 
by  the  simple  edge  and  the  previously  detected  edge  points, 
it  could  steer  around  the  breast  boundary  and  was  less  prone 
to  diversion  by  noise  and  artifacts.  The  accuracy  of  the 
boundary  tracking  technique  was  evaluated  in  our  previous 
study23  by  quantifying  the  root-mean-square  differences  be¬ 
tween  the  detected  and  manually  identified  breast  bound¬ 
aries.  In  the  current  study,  the  performance  of  the  boundary 
tracking  technique  for  this  data  set  was  determined  by  super¬ 
imposing  the  detected  boundary  on  the  breast  image  and  vi¬ 
sually  judged  if  the  detected  boundary  coincided  with  the 
perceived  breast  boundary.  The  breast  image  and  its  bound¬ 
ary  were  displayed  by  appropriately  adjusting  the  contrast 
and  brightness.  Incomplete,  jagged  and  mistracked  bound¬ 
aries  were  considered  incorrect  tracking. 

The  unexposed  film  area  around  the  film  edges  was  de¬ 
tected  automatically.  After  the  breast  boundary  was  found,  a 
region  growing  algorithm  was  used  to  fill  the  enclosed  breast 
region.  The  result  was  a  binary  map  that  distinguished  the 
breast  region  from  the  background  areas.  An  example  of  the 
tracked  breast  boundary  and  the  breast  binary  map  is  shown 
in  Figs.  l(a)-l(c). 

For  the  MLO  view  mammograms,  an  additional  step  has 
to  be  performed  for  segmentation  of  the  pectoral  muscle.  The 
initial  edge  in  the  pectoral  region  was  found  as  the  maximum 
gradient  point  by  a  line-by-line  gradient  analysis  from  the 
chest  wall  to  the  breast  boundary.  The  false  pectoral  muscle 
edge  points  were  discarded  by  an  edge  validation  process. 
First,  a  straight  line  was  fitted  to  the  initial  edge  points,  and 
the  points  that  did  not  lie  close  to  the  fitted  line  were  re¬ 
moved.  Second,  the  remaining  edge  points  that  were  con¬ 
nected  were  identified  by  an  8-connectivity  criterion.  An 
edge  segment  was  removed  if  its  direction  was  inconsistent 
with  the  pectoral  edge  direction  relative  to  the  breast  image. 
Finally,  a  second  order  curve  was  fitted  to  the  remaining 
edge  points  to  separate  the  pectoral  muscle  from  the  breast 
region.  The  pixels  in  the  pectoral  muscle  region  were  ex¬ 
cluded  from  the  histogram  analysis  and  breast  area  calcula¬ 
tion.  The  accuracy  of  the  pectoral  muscle  detection  was  also 
judged  visually  in  this  study,  similar  to  the  method  used  for 
the  breast  boundary  described  above.  Figure  1  shows  the 
pectoral  muscle  trimming  result  for  an  MLO  view  mammo¬ 
gram. 


To  facilitate  histogram  analysis,  a  dynamic  range  com¬ 
pression  method  was  developed  to  reduce  the  gray  level 
range  of  the  histograms.  With  our  digitization,  the  gray  lev¬ 
els  of  the  dense  tissue  are  higher  than  those  of  the  adipose 
tissue.  Because  of  variations  in  exposure  condition  and 
breast  thickness  near  the  periphery,  the  gray  level  distribu¬ 
tion  corresponding  to  the  breast  parenchymal  pattern  is  su¬ 
perimposed  on  a  low  frequency  background  that  mainly  rep¬ 
resents  the  global  variations  in  exposure.  This  low  frequency 
background  distorts  the  characteristic  features  of  the  histo¬ 
gram  due  to  the  density  pattern.  To  reduce  the  distortion,  an 
adaptive  dynamic  range  compression  technique  was  applied 
to  the  breast  image.  For  a  given  breast  image,  F(x,y),  which 
contains  low  frequency  background  and  higher  frequency 
breast  tissue  structures,  a  smoothed  image,  FB(x,y),  was 
obtained  by  applying  a  large-scale  box  filter  to  F(x,y)  to 
remove  the  high  frequency  components  while  retaining  the 
low  frequency  components.  The  image  FB(x,y)  was  then 
compressed  by  a  scale  factor  k : 

Fc(x,y)  =  kFB(x,y).  (1) 

To  reconstruct  the  high  frequency  components,  Fc(x,y), 
was  subtracted  from  a  constant  gray  level  G,  and  added  to 
the  original  image,  F(x,y): 

FD(x,y)  =  G~Fc(x,y),  (2) 

FE(x,y)=FD(x,y)+F(x,y).  (3) 

Histogram  analysis  was  applied  to  the  dynamic-range- 
compressed  image  FE(x,y).  Figure  2  shows  an  example  of 
the  resulting  images  and  gray  level  histograms  obtained  from 
this  procedure,  where  the  size  of  box  filter  is  35X35,  the 
scale  factor  k  is  0.5,  and  the  constant  gray  level  G  is  the 
maximum  gray  level  of  the  compressed  image  Fc(x>y).  The 
values  of  these  parameters  were  chosen  experimentally  as  a 
balance  between  reducing  the  dynamic  range  and  preserving 
the  image  features  in  the  compressed  image. 

C.  Breast  density  segmentation  and  estimation 

A  rule-based  threshold  technique  was  developed  to  seg¬ 
ment  the  dense  areas  from  the  breast  background.  The  histo¬ 
gram  of  the  breast  region  on  the  dynamic-range-compressed 
mammogram  was  generated  and  smoothed.  The  histograms 
of  these  images  in  the  database  were  analyzed  to  formulate 
an  automatic  thresholding  routine.  The  histograms  were 
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Fig.  2.  (a)  A  typical  mammogram  from  our  image  database;  (b)  the  low  frequency  image  Fa{ x,y)  obtained  by  an  35x  35  box  filter;  (c)  the  compressed  image 
(d)  the  inverted  image  FD(x,y);  (e)  the  enhanced  image  F£(*,y);  (0  the  gray  level  histogram  within  the  breast  region  of  the  original  image  F(jc,y); 
and  (g)  the  gray  level  histogram  of  the  breast  region  of  the  enhanced  image  FE{x,y). 


grouped  into  four  classes  based  on  the  characteristic  shapes 
of  their  histograms.  It  was  observed  that  the  grouping  corre¬ 
sponded  approximately  to  the  four  BI-RADS  breast  density 
ratings:  Class  I  corresponded  to  breasts  of  almost  entirely  fat. 
Class  II  corresponded  to  scattered  fibroglandular  densities, 
Class  ID  corresponded  to  heterogeneously  dense  and  Class 
IV  corresponded  to  extremely  dense  breasts.  Examples  of 
typical  histograms  for  these  four  classes  are  shown  in  Fig.  3. 
The  histograms  seemed  to  follow  two  basic  patterns.  In  one 
pattern,  there  was  only  one  dominant  peak,  which  repre¬ 
sented  most  of  the  breast  structures  in  the  breast  region.  In 
the  other  pattern,  in  addition  to  a  large  peak  in  the  histogram, 
there  was  one  or  two  smaller  peaks  on  the  right  or  left  side  of 
the  large  peak.  In  a  majority  of  the  cases,  the  smaller  peak 
was  distinguishable  from  the  large  one  when  the  random 
fluctuation  on  the  histogram  was  smoothed. 


1 .  Peak  detection  and  feature  description 

The  gray  level  histogram  within  the  breast  area  was  gen¬ 
erated  and  normalized,  and  passed  through  an  averaging  win¬ 
dow  to  smooth  out  the  random  fluctuations.  We  estimated 
the  window  size  to  be  in  the  range  of  30  to  50  gray  levels  by 
experimentally  evaluating  the  histogram  shapes  and  density 
segmentation  at  different  window  sizes.  Too  small  a  window 
size  cannot  smooth  out  the  fluctuation  and  too  large  a  win¬ 
dow  size  will  blur  the  useful  features.  A  window  size  of  30 
was  used  in  this  study.  The  second  derivative  of  every  point 
on  the  histogram  curve  was  computed.  An  example  of  the 
histogram  and  its  second  derivative  curve  are  shown  in  Fig. 
4.  The  zero  crossing  locations  were  detected  by  scanning  for 
the  positive-to-negative  and  negative- to-positive  changes  on 
the  latter  curve.  If  the  second  derivative  was  negative  be- 
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Fig.  3.  Four  typical  classes  of  histograms  and  the  setting  of  gray  level  interval  [gx  ,g2]  for  the  threshold  calculation. 


tween  two  zero  crossing  points,  it  indicated  that  a  peak  ex¬ 
isted  between  these  two  points  on  the  histogram.  Normally, 
as  shown  in  Fig.  4,  a  peak  included  the  peak  point  P0  and 
two  valley  points  Pj  and  P2  located  on  the  two  sides  of  the 
peak  point.  The  peak  point  P0  was  determined  by  searching 
for  the  maximum  histogram  value  between  the  zero  crossing 
points  Z2  and  Z3 ,  and  the  Px  and  P2  points  were  obtained  by 
searching  for  the  point  with  minimum  histogram  value  be¬ 
tween  zero  crossing  points  Z\  ,Z2  and  Z3  ,Z4,  respectively. 

The  following  peak  features  can  be  defined  by  peak  point 
P0  and  valley  points  P{  and  P2 : 

1  P2 

Energy:  E=  -  2  /(')*/('),  (4) 

A  i=p  j 

1  P° 

left-side  energy:  EL= j  X  /(*)*/(*).  (5) 

A  i=p 


Gray  level 

Fig.  4.  The  gray  level  histogram  (solid  curve)  and  the  second  derivative 
(dot)  curve.  P0  is  the  peak  point,  Pj  and  P2  are  the  valley  points  of  the  peak 
on  the  two  sides  of  the  peak  point  P0 .  Points  Zj ,  Z2 ,  Z3  and  Z4  are  zero 
crossing  points  on  the  second  derivative  curve,  which  are  used  for  searching 
the  points  P0 ,  Pj  and  P2. 


1  ^2 

right-side  energy:  ER= —  2  (6) 

A  i=Po 

likelihood:  L=EIE' ,  (7) 

where  /( ■ )  is  the  histogram,  A  is  the  total  energy  of  the 
entire  histogram  and  A  is  the  maximum 

gray  level  of  the  histogram.  Ef  is  the  energy  calculated  by 
approximating  the  histogram  in  the  interval  [Pi,P2]  using 
two  straight  lines,  P\Pq  P0P2  .  The  energy  E  of  the  peak 
is  used  to  compare  the  sizes  of  the  peaks  on  the  histogram, 
higher  energy  means  bigger  size  of  the  peak.  Eh  and  ER  split 
the  energy  E  into  two  parts  from  the  peak  point  for  calculat¬ 
ing  the  ratio  of  the  energy  in  these  two  parts.  The  likelihood 
L  describes  how  close  the  real  peak  is  to  the  triangle  repre¬ 
sented  by  the  three  points  P0,  Px  and  P2 . 

2.  Rule-based  histogram  classification 

A  rule-based  histogram  classifier  was  developed  to  clas¬ 
sify  the  gray  level  histogram  of  the  breast  area  into  four 
classes.  As  shown  in  Fig.  3,  a  typical  Class  I  breast  is  almost 
entirely  fat,  it  has  a  single  narrow  peak  on  the  histogram. 
Class  II  has  scattered  fibroglandular  densities,  it  has  two 
peaks,  other  than  the  tail  part  on  the  left,  on  the  histogram, 
with  the  smaller  peak  on  the  right  of  the  bigger  one.  Class  m 
is  heterogeneously  dense,  it  also  has  two  peaks,  but  the 
smaller  peak  is  on  the  left  of  the  bigger  one.  Class  IV  is 
extremely  dense,  which  has  a  single  dominant  peak  on  the 
histogram,  but  it  is  wider  compared  with  the  peak  in  the 
Class  I  histogram,  and  a  second  small  peak  sometimes  oc¬ 
curs  to  the  left  of  the  main  peak. 

The  classification  is  performed  in  two  steps.  In  the  first 
step,  the  computer  determines  whether  there  is  only  one 
single  peak  in  the  histogram.  The  biggest  peak  (main  peak) 
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PM  and  its  location  are  detected  by  comparing  the  energy  of 
the  peaks  on  the  histogram.  The  single  peak  feature  is  mainly 
determined  by  the  energy  E  under  the  main  peak  and  the 
features  EL  and  ER.  If  the  histogram  is  found  to  have  a 
single-peak  pattern,  in  general,  a  narrow  peak  corresponds  to 
very  fatty  breast  (Class  I),  and  a  wider  peak  corresponds  to 
very  dense  breast  (Class  IV).  However,  in  some  cases,  the 
histogram  of  these  two  classes  is  very  similar,  as  discussed 
below  (Fig.  9),  and  it  is  difficult  to  distinguish  them  by  their 
gray  level  histogram  distributions.  Two  additional  image  fea¬ 
tures  were  analyzed  to  classify  very  fatty  and  very  dense 
breasts.  One  feature  is  the  gray  level  standard  deviation  (Std) 
in  the  entire  breast  area,  defined  as 

{  1  \  1/2 

Std=  -  2  l(f(x,y)-/(x,y))2  ,  (8) 

\  Z’  x  e  MAP  y  e  MAP  J 

where  MAP  is  the  breast  binary  map  region,  N  is  the  pixel 
numbers  within  MAP.  Another  feature  is  the  number  of 
single  pixels  and  single  pixel-size  holes  (NSH)  counted  in 
the  breast  area  of  a  segmented  binary  image  using  the  big¬ 
gest  histogram  peak  point  PM  as  a  threshold.  For  a  very  fatty 
mammogram,  the  breast  mainly  consists  of  a  fatty  back¬ 
ground  with  some  fibrous  structures  and  fibroglandular  tissue 
scattered  in  the  breast  area.  The  NSH  value  was  found  to  be 
larger  (greater  than  50  pixels  on  average),  and  Std  smaller 
(less  than  500  on  average),  compared  with  a  mammogram  of 
a  very  dense  breast. 

In  the  second  step,  if  the  histogram  is  found  to  have  more 
than  one  peak,  decision  rules  are  used  to  decide  if  the  second 
major  peak  is  on  the  left  side  or  on  the  right  side  of  PM  by 
the  features  £,  ELy  ER  and  L,  and  the  relative  position  of  the 
two  peaks.  If  the  second  major  peak  is  on  the  right,  then  the 
histogram  is  classified  to  be  Class  II;  otherwise,  it  is  classi¬ 
fied  to  be  Class  HI. 

3 .  Gray  level  thresholding 

Gray  level  thresholding  is  essentially  a  pixel  classification 
problem.  Its  objective  is  to  classify  the  pixels  of  a  given 
image  into  two  classes:  one  includes  pixels  with  gray  values 
that  are  below  or  equal  to  a  certain  threshold;  the  other  in¬ 
cludes  those  with  gray  values  above  the  threshold.  Thresh¬ 
olding  is  a  popular  tool  for  image  segmentation,  a  variety  of 
techniques  have  been  proposed  over  the  years.  In  our  study, 
two  threshold  selection  methods  are  used:  one  is  the  Dis¬ 
criminant  Analysis  (DA)  method24  and  the  other  is  the  Maxi¬ 
mum  Entropy  Principle  (MEP)  based  method.25  The  DA 
method  assumes  that  the  image  gray  levels  can  be  classified 
into  two  classes  by  a  threshold.  To  estimate  the  threshold,  a 
discriminant  criterion  based  on  the  within-class  variance  and 
between-class  variance  is  introduced.  An  optimal  threshold 
is  selected  by  the  discriminant  criterion  to  maximize  the 
separability  of  the  resultant  classes  in  terms  of  gray  levels. 
This  method  is  well-suited  for  the  cases  where  the  gray  level 
histogram  is  bimodal.  In  an  ideal  situation,  the  histogram  has 
a  deep  and  sharp  valley  between  the  two  peaks  representing 
objects  and  background,  respectively,  and  the  optimum  cor¬ 
responds  to  the  gray  level  at  the  bottom  of  this  valley.  A 
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more  detailed  description  of  the  DA  method  can  be  found  in 
Appendix  A. 

For  the  MEP  method,  the  optimal  threshold  value  is  de¬ 
termined  by  maximizing  the  a  posteriori  entropy  subject  to 
certain  inequality  constraints  that  are  derived  by  means  of 
special  measures  characterizing  the  uniformity  and  the  shape 
of  the  regions  in  the  image.  As  is  well-known  26  the  maxi¬ 
mum  a  posteriori  probability  can  serve  as  a  criterion  to  se¬ 
lect  a  priori  probability  distributions  when  very  little  is 
known  about  the  probability  distribution.  Compared  with  the 
DA  method,  MEP  can  provide  a  better  thresholding  result  if 
the  gray  level  histogram  does  not  have  a  bimodal  distribu¬ 
tion.  A  more  detailed  description  of  the  MEP  method  can  be 
found  in  Appendix  B. 

The  gray  level  histograms  of  the  mammograms  in  our 
study  are  very  complex,  the  histogram  may  be  unimodal, 
bimodal  or  multi-modal.  It  is  difficult  to  select  an  appropriate 
threshold  by  one  general  threshold  selection  method.  There¬ 
fore,  we  combined  both  the  DA  and  the  MEP  methods,  to 
select  a  threshold  according  to  the  characteristic  features  of 
the  histogram  that  has  been  classified  into  one  of  the  four 
classes.  Suppose  /(g)  is  the  gray  level  histogram  of  the 
breast  area.  Let  7,=  Method(/(g)|g1<g<g2)  represent  the 
threshold,  T,  that  is  selected  by  use  of  Method  in  the  interval 
[gi  ,£2]  of  the  histogram  /(g),  where  Method  can  be  either 
the  DA  or  MEP  method.  The  settings  of  the  interval  [g]  ,g2] 
for  the  four  classes  are  discussed  below  and  shown  in  Fig.  3. 

Class  I:  The  histogram  is  unimodal  so  that  the  threshold  is 
selected  as 

r=MEP(/(g)|g1<5<g2), 

where,  gj  is  the  main  peak  point;  g2  is  the  valley  point  on 
the  right  side  of  main  peak. 

Class  II:  The  histogram  is  not  unimodal  and  the  histogram 
is  classified  as  Class  II;  the  threshold  is  selected  by  averaging 
two  thresholds  that  are  computed  in  two  different  intervals  of 
the  histogram  by  the  DA  method: 

71  =  DA(/(g)|g>g1), 

72  =  DA(/(g)|g>g2), 

T=(T}  +  T2)/2, 

where  gj  is  the  valley  on  the  left  of  the  main  peak;  g2  is  the 
main  peak  point. 

Class  HI:  The  histogram  is  not  unimodal;  there  are  two 
possibilities  in  the  histogram  distribution:  there  is  a  valley 
between  the  main  peak  and  its  left  side  peak,  as  shown  in 
Fig.  3,  or  no  obvious  valley  exists  between  the  main  peak 
and  its  left  side  peak.  In  two  different  intervals  of  the  histo¬ 
gram,  two  thresholds  are  computed  as 

7'i  =  DA(/(g)|g1<g<g2), 

T2  =  DMf(g)\g[<g<g2), 

where  gj  is  the  left  valley  point  of  the  left-side  peak  (PLM) 
of  the  main  peak,  g  [  is  the  peak  point  of  PLM  and  g2  is  right 
valley  point  of  the  main  peak.  If  there  is  an  obvious  valley, 
T=(T{  +  T2)I2,  otherwise  T=TX. 
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Class  IV:  Since  the  histogram  is  considered  unimodal,  the 
threshold  is  computed  by  the  MEP  method,  T 
=  MEP(f(g)\gl<g<g2),  where,  gx  is  the  left  valley  point 
of  the  main  peak;  g2  is  the  main  peak  point. 

D.  Radiologists’  segmentation  of  dense  breast  tissue 

In  order  to  evaluate  the  accuracy  of  the  computer  segmen¬ 
tation  method,  the  computer  segmentation  results  were  com¬ 
pared  to  radiologists’  manual  segmentation  in  the  data  set  of 
65  patient  cases.  Details  of  the  observer  study  for  estimation 
of  the  breast  density  and  statistical  analysis  of  the  results 
were  discussed  elsewhere.27  Briefly,  a  graphical  interface 
was  developed  for  displaying  the  mammograms  and  record¬ 
ing  the  observer’s  evaluation.  The  CC-view  and  MLO-view 
mammograms  for  a  given  breast  were  displayed  side-by- 
side;  a  radiologist  observer  examined  the  mammograms  and 
gave  a  BI-RADS  rating  and  a  visual  estimation  of  the  per¬ 
cent  breast  density  with  10%  increments.  After  the  subjective 
evaluation,  each  view  was  displayed  sequentially,  together 
with  the  histogram  of  the  dynamic-range-compressed  image. 
The  radiologist  would  interactively  choose  a  threshold  by 
moving  a  slider  along  the  abscissas  of  the  histogram  plot. 
The  segmented  binary  image,  displayed  side-by-side  with  the 
mammogram,  would  change  instantaneously  when  the 
threshold  was  changed.  The  radiologist  could  inspect  if  the 
segmented  area  corresponded  to  the  dense  area  on  the  mam¬ 
mogram.  Once  the  radiologist  was  satisfied  with  the  segmen¬ 
tation  of  the  dense  area,  the  gray  level  threshold  and  the 
percent  dense  area  derived  from  this  threshold  were  re¬ 
corded.  The  display  then  moved  to  the  next  view  of  the  same 
breast  for  evaluation.  The  mammograms  of  the  other  breast 
for  the  same  patient  would  then  be  displayed  and  evaluated 
in  the  same  way.  The  entire  process  was  repeated  for  each 
patient  until  all  patients  in  the  data  set  were  evaluated. 

Five  MQSA-approved  radiologists  participated  in  the  ex¬ 
periment.  To  familiarize  the  radiologists  with  the  procedures 
and  to  assist  them  in  their  visual  estimation  of  the  percent 
breast  density,  we  had  them  trained  on  a  separate  set  of  25 
patient  cases  prior  to  the  evaluation  of  the  actual  data  set. 
During  the  training  session,  the  computer  displayed  the  per¬ 
cent  breast  dense  area  to  the  radiologist,  which  was  obtained 
by  the  radiologist’s  interactive  thresholding  of  the  image. 
The  radiologist  could  then  compare  the  manually  segmented 
percentage  with  their  visually  assessed  percent  density  for 
the  image.  This  feedback  helped  “calibrate”  the  radiolo¬ 
gists’  visual  estimates  of  the  percent  dense  breast  area.  The 
percent  dense  area  obtained  by  interactive  thresholding  was 
not  displayed  during  the  actual  study. 

III.  RESULTS 

An  example  of  a  typical  mammogram  from  each  of  the 
four  classes  and  its  corresponding  enhanced  image,  its  histo¬ 
gram,  the  selected  threshold  and  the  segmented  image  are 
shown  in  Figs.  5 (a) -5(d),  respectively. 

The  average  percent  breast  density  obtained  from  manual 
segmentation  by  the  five  trained  radiologists  for  each  mam¬ 
mogram  was  used  as  the  “true  standard”  of  the  percent 
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breast  density  for  that  mammogram.  The  breast  region  was 
segmented  by  the  breast  boundary  tracking  technique,  and 
the  pectoral  muscle  was  trimmed  for  the  MLO-view  mam¬ 
mograms.  The  breast  boundary  was  accurately  tracked  on 
92.3%  (240/260)  of  the  mammograms,  and  the  pectoral 
muscle  was  correctly  trimmed  on  74.6%  (97/130)  of  the 
MLO  views.  The  histograms  of  6%  (8  CC  views  and  8  MLO 
views)  of  the  breast  regions  did  not  exhibit  the  typical  char¬ 
acteristic  features  of  the  four  classes  and  were  misclassified 
by  the  computer,  resulting  in  poor  segmentation  of  the  dense 
region. 

Figure  6  shows  a  comparison  of  the  percent  breast  density 
visually  estimated  by  radiologists  against  the  true  standard 
for  the  94%  of  the  260  mammograms  that  were  classified 
correctly  by  the  computer.  Table  I  summarizes  the  compari¬ 
son  of  the  radiologists’  visual  estimates  with  the  true  stan¬ 
dard.  The  “difference”  between  the  estimated  %  breast  den¬ 
sity  and  the  true  standard  was  calculated  for  each  case,  and 
the  mean  and  the  standard  deviation  of  this  difference  over 
all  cases  were  estimated  for  each  radiologist  and  shown  in 
the  table.  Therefore,  the  mean  difference  was  the  average 
bias  of  the  estimated  %  breast  density  from  the  true  standard 
over  all  images  in  the  data  set.  It  can  be  seen  that  almost  all 
radiologists  had  a  positive  bias,  on  average,  when  they  visu¬ 
ally  estimated  mammographic  density,  except  for  Radiolo¬ 
gist  5  who  had  a  small  negative  average  bias  on  the  CC-view 
reading.  For  a  given  radiologist,  the  over-estimation  in¬ 
creased  as  the  breast  density  increased.  Although  the  corre¬ 
lation  coefficients  were  high,  ranging  from  0.90  to  0.95,  the 
deviations  from  the  diagonal  line  were  systematic.  The  aver¬ 
age  bias  from  the  true  standard  varied  from  less  than  1%  to 
11%,  depending  on  the  radiologist.  The  root- mean-square 
(RMS)  errors  of  the  five  radiologists  relative  to  the  true  stan¬ 
dard  ranged  from  7.5%  to  16.3%. 

Figure  7  shows  the  comparison  of  the  percent  breast  den¬ 
sity  between  the  computer  segmentation  and  the  true  stan¬ 
dard  for  the  94%  of  mammograms  whose  histograms  were 
considered  to  be  correctly  classified.  There  was  a  trend  of 
over-estimation  in  the  very  fatty  breasts.  In  the  medium 
dense  range,  the  variances  from  the  true  standard  were  high. 
Some  images  had  a  large  deviation  from  the  diagonal  line, 
indicating  that  the  threshold  was  incorrectly  determined. 
Table  II  summarizes  the  comparison  between  the  computer 
performance  and  the  true  standard.  For  the  CC  views  with 
correct  histogram  classification,  the  correlation  between  the 
computer-estimated  percent  dense  area  and  the  true  percent 
breast  density  was  0.94,  and  between  the  computer  and  the 
radiologists’  average  visual  estimate  was  0.87  (not  plotted). 
These  correlation  coefficients  were  0.91  and  0.82,  respec¬ 
tively,  for  the  MLO  views  with  correct  classification.  Al¬ 
though  the  correlation  coefficients  of  the  computer  segmen¬ 
tation  with  the  true  standard  were  not  better  than  those  of  the 
visual  estimates,  the  average  biases  of  the  computer  segmen¬ 
tation  from  the  true  standard  were  less  than  2%,  which  were 
substantially  less  than  those  of  visual  estimates  (Table  I). 
This  indicates  that  computerized  segmentation  is  a  good  al¬ 
ternative  to  manual  segmentation  although  variances  of  the 
automated  method  will  need  to  be  further  reduced.  The  RMS 
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Fig.  5.  Four  classes  of  typical  mammograms  and  corresponding  enhanced  and  segmented  image,  histogram  and  threshold. 


errors  of  the  computer  segmentation  were  also  less  than 
those  of  the  radiologists’  visual  estimates,  at  6.1%  and  7.2%, 
respectively,  for  the  CC  view  and  MLO  view,  when  the  his¬ 
tograms  were  correctly  classified.  The  biases  and  RMS  errors 
for  the  different  subsets  of  images  are  also  shown  in  Table  EL 
It  can  be  seen  that  correct  histogram  classification  was  the 
most  important  factor  in  reducing  the  biases  and  the  RMS 
errors.  The  contributions  by  breast  boundary  detection  and 
pectoral  muscle  segmentation  were  minor,  on  average,  for 
improving  the  estimation  of  the  percent  dense  breast  area. 

Figure  8  shows  the  comparison  of  the  individual  radiolo¬ 
gists’  manual  segmentation  against  the  true  standard.  For  CC 
views,  the  RMS  difference  in  the  percent  breast  density  be¬ 


tween  an  individual  radiologist’s  manual  segmentation  and 
the  true  standard  varied  from  2.9%  to  5.9%  among  the  five 
radiologists.  For  MLO  views,  the  RMS  difference  varied 
from  2.8%  to  6.2%.  The  average  biases  of  the  five  radiolo¬ 
gists  ranged  from  -2.8%  to  2.2%  for  the  CC  views  and  from 
-3.1%  to  3.0%  for  the  MLO  views.  The  maximum  biases  of 
the  five  radiologists  varied  from  4.4%  to  22.6%  for  the  CC 
views  and  from  5.2%  to  23%  for  the  MLO  views. 

The  five  radiologists  provided  BI-RADS  density  ratings 
for  each  breast.  Although  the  BI-RADS  ratings  exhibited 
large  inter-observer  variations,20  it  is  interesting  to  compare 
the  computer’s  histogram  classification  with  the  BI-RADS 
ratings.  Since  there  were  260  images,  each  with  5  radiolo- 
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MLO  view 


Fig.  6.  A  comparison  of  the  percent  breast  density  between  five  radiologists’ 
visual  estimates  and  the  true  standard.  The  dashed  line  represents  the  linear 
regression  of  all  data  points  on  the  plot.  The  MLO  view  is  shown.  The  trend 
for  the  CC  view  is  similar. 

gists’  ratings,  there  were  a  total  of  1300  rating  comparisons. 
The  comparison  of  the  computer  and  the  radiologists’ 
BI-RADS  ratings  is  shown  in  Table  IE.  It  was  found  that 
87.4%of  Class  I  classification  have  BI-RADS  ratings  1  or  2, 
92.0%  of  Class  II  classifications  have  density  ratings  2  or  3, 
83.4%  of  Class  HI  classifications  have  density  ratings  3  or  4 


and  57.1%  of  Class  IV  classifications  have  density  rating  4. 
More  detailed  analysis  of  the  variability  of  radiologists’  BI¬ 
RADS  ratings  was  discussed  by  Martin  et  ol. 21 

IV.  DISCUSSION 

Radiologists  routinely  estimate  mammographic  breast 
density  using  the  four  BI-RADS  categories.  In  studies  that 
require  breast  density  estimation,  radiologists’  visual  esti¬ 
mates  of  mammographic  density  were  often  used  as  the  den¬ 
sity  measure.  Our  observer  study  indicates  that  interobserver 
variation  between  the  BI-RADS  ratings  of  five  experienced 
radiologists  ranged  from  —1  to  +1.  The  subjectively  esti¬ 
mated  percent  dense  area  can  deviate  from  the  true  standard 
by  as  much  as  40%,  as  shown  in  Fig.  6.  These  results  indi¬ 
cate  the  need  to  develop  an  objective  method  for  the  estima¬ 
tion  of  mammographic  breast  density  in  order  to  improve  the 
accuracy  and  reproducibility  of  the  estimation.  A  computer¬ 
ized  image  analysis  method  for  mammographic  breast  den¬ 
sity  estimation  will  be  a  useful  tool  for  study  of  breast  cancer 
risk  factors  and  for  monitoring  the  change  of  breast  cancer 
risk  with  preventive  or  interventional  treatments. 

In  this  study,  we  used  the  average  of  the  percent  breast 
area  obtained  with  interactive  thresholding  by  five  experi¬ 
enced  radiologists  as  the  true  standard.  The  gray  level  thresh¬ 
olding  method  used  in  this  study  could  achieve  a  reasonable 
segmentation  of  the  dense  areas  on  the  mammogram  because 
the  image  was  preprocessed  with  dynamic  range  compres¬ 
sion.  The  image-based  analysis  of  breast  density  will  not 
provide  the  actual  percentage  of  fibroglandular  tissue  in  the 
breast  volume.  However,  the  previous  studies  that  estab- 


Table  I.  A  comparison  of  the  radiologists’  visual  estimate  of  mammographic  breast  density  with  the  true 
standard.  The  “difference”  was  defined  as  the  difference  between  the  estimated  %  breast  density  and  the  true 
standard  for  each  case,  and  the  mean  and  the  standard  deviation  of  this  difference  are  tabulated. 


Image  subsets 

No.  of 
images 

Radiologist 

Correlation 

RMS 

error 

Mean 

difference 

Std.  dev.  of 
difference 

CC  view: 

All 

130 

Rad.  1 

0.942 

13.3% 

6.9% 

11.5% 

Rad.  2 

0.931 

14.5% 

9.8% 

10.7% 

Rad.  3 

0.923 

13.3% 

6.3% 

11.8% 

Rad.  4 

0.934 

7.5% 

2.9% 

7.0% 

Rad.  5 

0.901 

9.6% 

-1.4% 

9.6% 

Histogram 

122 

Rad.  1 

0.946 

13.7% 

7.2% 

11.3% 

correctly 

Rad.  2 

0.936 

14.7% 

10.3% 

10.8% 

classified 

Rad.  3 

0.929 

14.2% 

6.7% 

11.6% 

Rad.  4 

0.929 

7.7% 

3.1% 

7.1% 

Rad.  5 

0.900 

9.7% 

-1.3% 

9.4% 

MLO  view: 

All 

130 

Rad.  1 

0.933 

14.5% 

8.3% 

12.0% 

Rad.  2 

0.914 

16.1% 

11.2% 

11.5% 

Rad.  3 

0.915 

14.4% 

7.7% 

12.2% 

Rad.  4 

0.919 

8.8% 

4.3% 

7.7% 

Rad.  5 

0.910 

9.2% 

0.1% 

9.2% 

Histogram 

122 

Rad.  1 

0.932 

15.0% 

8.3% 

12.0% 

correctly 

Rad.  2 

0.914 

16.3% 

10.9% 

11.4% 

classified 

Rad.  3 

0.919 

14.7% 

7.8% 

12.2% 

Rad.  4 

0.916 

9.0% 

4.3% 

7.7% 

Rad.  5 

0.909 

9.4% 

0.3% 

9.2% 
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Fig.  7.  A  comparison  of  the  percent  breast  density  between  the  computer 
segmentation  and  the  true  standard.  The  dashed  line  represents  the  linear 
regression  of  the  data  on  the  plot,  (a)  CC  view,  (b)  MLO  view. 


lished  the  correlation  between  breast  density  and  breast  can¬ 
cer  risk  were  all  based  on  mammographic  density.  This  in¬ 
dicated  that  mammographic  density  is  a  sufficiently  sensitive 
marker  for  breast  cancer  risk,  although  it  may  be  less  accu¬ 
rate  than  volumetric  density.  An  actual  measurement  of  the 
percentage  of  fibroglandular  tissue  volume  in  the  breast,  for 
example,  by  x-ray  penetration  with  correction  for  scatter  and 
beam  hardening,  is  difficult  because  it  requires  accurate 
x-ray  sensitometry  or  phantom  calibration  for  each  image. 
These  requirements  will  limit  its  use  to  a  few  laboratories 
that  have  specialized  equipment  and  expert  physicists.  Mag¬ 
netic  resonance  breast  imaging  can  also  provide  volume 
measurement  of  dense  tissue  but  it  is  expensive  and  not  eas¬ 
ily  accessible.  It  can  be  expected  that  the  estimation  of  mam- 
mographic  breast  density  by  a  computerized  image  analysis 
method  will  be  a  more  practical  and  viable  approach,  espe- 
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dally  when  direct  digital  mammography  becomes  more 
widely  used  in  the  future. 

Our  preliminary  study  indicates  that  breast  density  esti¬ 
mation  can  be  performed  automatically  and  accurately  (Fig. 
7).  Although  the  accuracy  of  our  current  algorithm  still  needs 
to  be  improved,  it  can  be  seen  that  the  computer  segmenta¬ 
tion  can  provide  an  estimate  of  the  percent  breast  density 
with  a  very  small  bias  (Table  II).  More  importantly,  com¬ 
puter  segmentation  will  be  more  reproducible  and  consistent 
than  visual  estimates.  This  will  improve  the  sensitivity  of 
studies  that  depend  on  evaluation  of  the  change  in  mammo¬ 
graphic  density  over  time  or  before  and  after  a  certain  treat¬ 
ment. 

In  this  study,  we  reduced  the  spatial  resolution  to  a  pixel 
size  of  800/i,mX800/tm  for  image  processing.  The  small 
matrix  size  of  the  reduced  images  improves  the  computa¬ 
tional  efficiency.  The  reduction  in  resolution  has  two  major 
effects:  reducing  the  image  noise  and  blurring  the  details. 
Since  the  significant  dense  tissue  in  the  breast  that  contrib¬ 
utes  to  the  parenchyma  is  relatively  large  compared  to  800 
fi m,  it  is  not  expected  that  processing  at  this  pixel  size  will 
have  a  strong  effect  on  the  accuracy  of  the  estimated  percent 
breast  density.  Differences  in  the  segmented  area  may  occur 
mainly  along  the  boundary  of  the  dense  tissue  region,  but  the 
effect  may  be  averaged  out  statistically  along  boundaries  of 
reasonable  lengths.  The  residual  errors  in  the  estimation  of 
the  dense  area  should  not  be  substantial  in  comparison  with 
the  inter-  and  intra-radiologists*  variations  in  their  manual 
segmentation. 

Successful  segmentation  of  dense  tissue  depends  strongly 
on  whether  a  mammogram  can  be  classified  correctly  into  a 
proper  class.  A  successful  classification  will  likely  result  in 
the  selection  of  a  near  optimal  threshold.  Conversely,  if  a 
mammogram  is  classified  into  a  wrong  class,  the  threshold 
will  be  selected  incorrectly.  For  the  mammograms  of  very 
fatty  breasts,  the  gray  level  histogram  has  the  characteristics 
of  Class  I,  which  contains  one  large  single  peak.  These  his¬ 
tograms  can  be  distinguished  relatively  easily  from  most  of 
the  other  classes  of  histograms  if  those  histograms  exhibit 
the  typical  features.  For  mammograms  of  BI-RADS  category 
2  or  3,  there  are  scattered  fibroglandular  or  heterogeneous 
densities  in  the  breast.  A  small  peak  may  be  located  on  the 
left  or  on  the  right,  or  on  both  sides  of  the  main  peak  on  the 
histogram.  The  histogram  could  be  classified  into  Class  I  if 
the  small  peak  is  not  large  enough  and  is  not  detected  as  a 
second  peak.  Otherwise,  it  would  be  classified  into  Class  II 
or  Class  HI,  depending  on  the  location  of  that  small  peak 
relative  to  the  main  peak  of  the  histogram.  For  the  two-peak 
pattern  histogram,  the  DA  threshold  selection  method  is  ro¬ 
bust  if  there  is  an  obvious  valley  between  the  two  peaks.  If 
the  valley  is  fiat  or  not  obvious,  averaging  the  two  thresholds 
obtained  by  the  DA  method  in  two  different  intervals,  as 
designed  for  this  study,  can  reduce  the  chance  of  calculating 
an  incorrect  threshold  that  differs  greatly  from  the  optimum, 
but  it  also  reduces  the  chance  of  finding  the  optimal  thresh¬ 
old.  Overall,  the  rules  designed  for  classification  of  the  two- 
peak  patterns  seem  to  perform  consistently  well  for  this  data 
set.  One  of  the  difficult  situations  is  to  distinguish  between 
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Table  II.  A  comparison  of  computer  segmentation  with  the  true  standard.  The  “difference”  was  defined  as  the 
difference  between  the  estimated  %  breast  density  and  the  true  standard  for  each  case,  and  the  mean  and  the 
standard  deviation  of  this  difference  are  tabulated. 


Image  subsets 

No.  of 
images 

Correlation 

RMS 

error 

Mean 

difference 

Std.  dev.  of 
difference 

CC  view: 

All 

130 

0.746 

12.3% 

1.3% 

12.3% 

Boundary  correctly  tracked 

120 

0.780 

11.4% 

1.4% 

11.4% 

Histogram  correctly  classified 

122 

0.943 

6.1% 

0.2% 

6.2% 

Boundary  and  histogram  correctly  done 

113 

0.953 

5.6% 

0.8% 

5.6% 

MLO  view: 

All 

130 

0.780 

11.6% 

1.9% 

11.5% 

Boundary  correctly  tracked 

120 

0.766 

11.9% 

2.1% 

11.7% 

Histogram  correctly  classified 

122 

0.914 

7.2% 

1.5% 

7.1% 

Pectoral  muscle  correctly  trimmed 

97 

0.733 

11.6% 

1.6% 

11.6% 

Boundary  and  histogram  correctly  done 

112 

0.912 

7.2% 

1.7% 

7.1% 

Boundary,  histogram  and  pectoral 

83 

0.891 

7.1% 

1.9% 

6.8% 

muscle  correctly  done 


Class  I  and  Class  IV,  when  the  histogram  of  a  very  dense 
breast  mimics  that  of  a  very  fatty  breast,  as  shown  in  Fig.  9. 
This  image  was  correctly  classified  with  the  additional  fea¬ 
tures,  Std  and  NSH.  However,  there  were  other  cases  that 
failed  in  spite  of  the  additional  criteria.  The  large  difference 
in  the  optimal  threshold  locations  between  these  two  classes 
will  lead  to  a  large  error  in  the  estimated  percent  breast  den¬ 
sity  if  the  histogram  is  misclassified.  Further  study  is  needed 
to  more  accurately  distinguish  these  two  classes. 

The  dynamic  range  reduction  technique  reduces  the  vari¬ 
ability  of  the  gray  level  histograms  and  enhances  their  char¬ 
acteristics.  This  pre-processing  facilitates  the  classification  of 
the  image  into  the  correct  class.  There  are  many  image 
smoothing  techniques  published  in  the  literature.  Low-pass 
filtering  with  a  box  filter  is  the  simplest  choice.  The  effec- 


MLO  view 


Fig.  8.  A  comparison  of  the  percent  breast  density  obtained  from  the  five 
radiologists’  manual  segmentation  with  their  average  for  the  same  mammo¬ 
grams.  The  MLO  view  is  shown.  The  trend  for  the  CC  view  is  similar. 


tiveness  of  background  correction  with  a  box  filtered  image 
depends  on  the  box  size.  We  found  that  a  35  X  35-pixel  filter 
is  a  good  balance  between  computation  time  and  the  capa¬ 
bility  to  remove  the  high  frequency  components.  The  sub¬ 
traction  of  the  low-pass  filtered  image  from  the  original  im¬ 
age  is  a  form  of  unsharp  masking.  The  breast  boundary  is 
generally  enhanced  as  shown  in  Fig.  2(e).  The  pixels  at  the 
enhanced  breast  boundary  contribute  a  small  peak  to  the  left 
tail  of  the  gray  level  histogram  of  the  breast  area.  Moreover, 
if  dense  tissue  is  present  close  to  the  breast  boundary,  it  may 
not  be  segmented  correctly  due  to  intensity  reduction.  Other 
low  frequency  estimation  techniques  such  as  wavelet  decom¬ 
position  will  be  investigated  in  future  studies. 

In  this  feasibility  study,  we  used  a  small  data  set  of  mam¬ 
mograms  to  develop  a  rule-based  classifier  for  the  histogram 
analysis.  Although  a  large  fraction  of  the  histograms  mani¬ 
fest  characteristic  features  that  can  be  grouped  into  four 
classes,  corresponding  approximately  to  the  four  BI-RADS 
breast  density  ratings,  there  are  many  exceptions.  One  such 
example  is  shown  in  Fig.  9.  This  causes  misclassification  and 
incorrect  thresholding  by  the  histogram  classifier.  It  will  be 


Table  HI.  A  comparison  of  computer  classification  and  radiologists’ 
BI-RADS  breast  density  ratings. 


Computer 

classification 

BI-RADS 

1 

BI-RADS 

2 

BI-RADS 

3 

BI-RADS 

4 

Total 

Class  I 

210 

262 

52 

16 

540 

(16.2%) 

(20.2%) 

(4%) 

(1.2%) 

(41.5%) 

Class  II 

0 

92 

184 

24 

300 

(0%) 

(7.1%) 

(14.2%) 

(1.8%) 

(23.1%) 

Class  in 

1 

52 

167 

100 

320 

(0.1%) 

(4%) 

(12.8%) 

(7.7%) 

(24.6%) 

Class  IV 

5 

12 

43 

80 

140 

(0.4%) 

(0.9%) 

(3.3%) 

(6.2%) 

(10.8%) 

Total 

216 

418 

446 

220 

1300 

(16.6%) 

(32.2%) 

(34.3%) 

(16.9%) 

(100%) 
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Fig.  9.  The  gray  level  histograms  of  two  mammograms  classified  by  radiologists  as  BI-RADS  rating  1  (upper  mammogram)  and  BI-RADS  rating  4  (lower 
mammogram).  The  shapes  of  the  histograms  are  very  similar  and  cannot  be  distinguished  by  our  current  histogram  analysis  method.  These  two  examples  were 
correctly  classified  with  the  additional  Std  and  NSH  criteria. 


necessary  to  investigate  if  other  classification  strategies  can 
be  more  effective  than  a  rule-based  method.  Furthermore,  we 
have  not  performed  a  systematic  study  to  optimize  the  many 
parameters  used  in  the  segmentation  algorithm.  Further  work 
will  be  required  to  investigate  the  dependence  of  the  segmen¬ 
tation  accuracy  on  the  various  parameters.  The  parameter 
selection  and  the  performance  of  the  computer  classifier  will 
have  to  be  improved  by  training  with  a  larger  data  set  and  its 
generalizability  evaluated  with  unknown  cases.  The  generali¬ 
zation  of  the  algorithm  to  images  acquired  with  other  digi¬ 
tizers  or  direct  digital  mammography  systems  will  also  need 
to  be  investigated. 

V.  CONCLUSION 

We  are  developing  an  image  analysis  method  for  auto¬ 
mated  segmentation  of  the  dense  area  from  mammograms 
and  estimation  of  the  percent  mammographic  density.  Our 
preliminary  study  indicates  the  feasibility  of  our  approach. 
The  computer-estimated  mammographic  breast  density  cor¬ 
relate  closely  with  the  average  manual  segmentation  by  five 
experienced  radiologists  and  the  average  bias  is  much  less 
than  that  of  the  radiologists’  visual  estimation.  We  have 
found  that  correct  classification  of  the  histogram  shapes  is 
the  most  crucial  step  in  our  approach.  The  histograms  of 
many  mammograms  have  distinctive  characteristics  that  can 
be  recognized  by  a  rule-based  classifier.  However,  some  his¬ 
tograms  deviate  from  these  rules  and  this  can  lead  to  mis- 
classifi cation.  A  further  investigation  will  be  needed  to  de- 
sign  more  robust  rules  or  classifiers  to  improve  the 
classification  accuracy.  Despite  these  limitations,  we  have 


demonstrated  in  this  preliminary  study  that  the  estimation  of 
mammographic  density  can  be  performed  efficiently  and  ac¬ 
curately  by  the  automated  image  analysis  tool.  The  fully  au¬ 
tomated  algorithm  can  provide  an  objective  and  reproducible 
quantitative  estimation  of  mammographic  breast  density  that 
is  expected  to  be  superior  to  subjective  visual  assessment  and 
comparable  to  manual  segmentation  by  radiologists. 
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APPENDIX  A:  GRAY-LEVEL 

THRESHOLDING— DISCRIMINANT  ANALYSIS  (DA) 

METHOD 

Suppose  the  probability  of  the  gray  level  n,  in  an  image 
with  L  gray  levels  can  be  estimated  as 

L 

Pi-tiilN,  N=2  (Al) 

i=  1 

If  the  pixels  in  the  image  are  classified  into  two  classes  C0 
and  Cx  by  the  threshold  k,  then  the  probabilities  of  class 
occurrence  and  the  class  mean  levels  are  given  by 
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k  L 

"0=2  Pi  =  P(k),  &>!=  2  Pi=l-P(k), 

1=1  i=k  + 1 

k 

/io =  2  iPil<o0=fi(k)/(o0, 

i=  1 

'‘■“A,  ",'/"'=T^w- 

where 


(A2) 


(A3) 


k  k  L 

^W=2  Pi.  P-(*)  =  2  iPi  and  /ir=2  ipt,  (A4) 

i=l  i=  1  i=l 


are  the  zeroth-  and  the  first-order  cumulative  moments  of  the 
histogram  up  to  the  kth  level,  and  the  total  mean  level  of 
original  image,  respectively. 

The  between-class  variance  is  defined  as 


o-2B(k)  =  o>0(fj.0- fjiT)2+  fiT)2 

/  x2  _\-P-TP(k)-fi(k)]2 

— o-ito-#*)  -  P(jfc)[i-P(I)T> 

and  the  optimal  threshold  A:*  is  given  by 

max  VbW- 

l^k^L 


(A5) 

(A6) 


APPENDIX  B:  GRAY-LEVEL 
THRESHOLDING— MAXIMUM  ENTROPY 
PRINCIPLE  (MEP)  METHOD 

Suppose  the  probability  of  the  gray  level  n(  in  an  image 
with  L  gray  levels  can  be  estimated  as 

L 

Pi=ni/N,  N =2  n,.  (Bl) 

1=1 

After  thresholding  the  image  by  threshold  k ,  the  a  posteriori 
probability  of  the  pixels  with  gray  level  value  less  than  k ,  is 
given  by 

k 

F(k)  =  J,Pi.  (B2) 

i=0 

And  the  a  posteriori  probability  of  all  those  pixels  with  val¬ 
ues  greater  than  or  equal  to  k  is  1  -F(k).  Thus  the  Shannon 
entropy  of  the  thresholded  image  is 

H(F(k))=  -F(k)log  F(k)-(\  -F(k)) log(l  -F(*)). 

(B3) 

The  optimal  threshold  k  maximizes  H(F(k)). 

^Correspondence:  Heang-Ping  Chan,  Ph.D.,  University  of  Michigan,  De¬ 
partment  of  Radiology,  1500.  Medical  Center  Drive,  UHB1F510B,  Ann 
Arbor,  Michigan  48109-0030.  Electronic  mail:  chanhp@umich.edu 
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ABSTRACT 

Studies  have  demonstrated  a  strong  correlation  between  mammographic 
breast  density  and  breast  cancer  risk.  Mammographic  breast  density  may 
therefore  be  used  as  a  surrogate  marker  for  monitoring  the  response  to 
treatment  in  studies  of  breast  cancer  prevention  or  intervention  methods. 

in  this  study,  we  evaluated  the  accuracy  of  using  mammograms  for 
estimating  breast  density  by  analyzing  the  correlation  between  the  percent 
mammographic  dense  area  and  the  percent  fibroglandular  tissue  volume 
obtained  from  MR  images.  A  data  set  of  37  patients  who  had  corresponding 
MR  images  and  mammograms  was  collected.  The  fibroglandular  tissue 
regions  in  the  MR  slices  were  segmented  by  a  semi-automatic  method  and 
the  percent  fibroglandular  tissue  volume  calculated.  Mammographic  breast 
density  was  estimated  by  an  automated  image  analysis  program- 
Mammographic  Density  ESTimator  (MDEST)  that  we  developed  previously.  It 
was  found  that  the  correlation  between  the  percent  dense  area  of  the  CC  and 
MLO  views  and  the  percent  fibroglandular  tissue  volume  on  MR  images  was 
0.93  and  0.91,  respectively,  with  a  mean  bias  of  4.4%.  The  high  correlation 
indicates  the  usefulness  of  mammographic  density  as  a  surrogate  for  breast 
density  estimation. 


MATERIALS  AND  METHODS 

DATA  SET 

37  patients,  corresponding  MR  breast  images  and  mammograms: 

MR  series:  Coronal  3DT1-weighted  pre-contrast  images 
3D  Spoiled  Gradient-Recalled  Echo  (SPGR) 

TE=3.3  ms;TR=10  ms,  Flip=40  deg,  matrix=256x128, 
FOV=28-32  cm  right/left,  14-16  cm  superior/inferior, 
scan  time=2  min  38  sec 
Coronal  sections  2-5  mm  thick,  32  slices; 

Full  volume  coverage  of  both  breasts  with  contiguous  slices 

Mammograms:  GE  DMR  mammography  system  (Mo/Mo,  Mo/Rh,  Rh/Rh) 
CC  and  MLO  views 


MR  BREAST  IMAGE  SEGMENTATION 


MAMMOGRAM  SEGMENTATION 


Breast  Boundary  Detection 


jjj  I  Background  Correction 


c  I  Gray  Level  Thresholding 


Morphological  Erosion 
for  exclusion  of  skin 


Fibroglandular  Tissue 
Volume  over  all  slices 


Percent  Fibroglandular 
Tissue  Volume 


Example  of  MR  images  from  one 
breast.  Although  the  MR  pulse 
sequence  Is  not  the  optimal  one  for 
separating  water  and  fat,  the 
fibroglandular  parenchyma  (-water) 
and  fatty  tissue  are  well  separated 
with  this  heavily  T1  -weighted 
acquisition.  Breast  boundaries  were 
automatically  detected  and  then 
manually  corrected  if  necessary. 
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Digitized  Mammogram 


Breast  Boundary  Detection 


Dynamic  Range 
Compression 


Histogram  Analysis 


Breast  Density 
Classification 


Automatic 

Gray  Level  Thresholding 


Percent  Dense  Area 


An  example  of  mammographic  density 
segmentation.  Two  views  (CC  and  MLO)  for 
each  breast  were  analyzed. 


CONCLUSION 

Although  the  %  projected  area  and  the  %  volume  of  fibroglandular  tissue  in  the 
breast  are  not  expected  to  be  linearly  related,  the  correlation  between  the  two 
estimates  was  found  to  be  very  high.  This  result  indicates  the  usefulness  of 
mammographic  density  as  a  surrogate  of  %  fibroglandular  tissue  volume  for  study  of 
the  correlation  between  breast  density  and  breast  cancer  risk.  The  computerized 
image  analysis  tool  (MDEST)  Is  useful  for  estimation  of  mammographic  density.  The 
automated  analysis  is  expected  to  contribute  to  the  understanding  of  the  relationship 
of  mammographic  density  to  breast  cancer  risk,  detection,  and  prognosis,  and  to  the 
prevention  and  treatment  of  breast  cancer. 
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Abstract.  Studies  have  demonstrated  a  strong  correlation  between 
mammographic  breast  density  and  breast  cancer  risk.  Mammographic  breast 
density  may  therefore  be  used  as  a  surrogate  marker  for  monitoring  the  response 
to  treatment  in  studies  of  breast  cancer  prevention  or  intervention  methods.  In 
this  study,  we  evaluated  the  accuracy  of  using  mammograms  for  estimating 
breast  density  by  analyzing  the  correlation  between  the  percent  mammographic 
dense  area  and  the  percent  glandular  tissue  volume  as  estimated  from  MR 
images.  A  data  set  of  37  patients  who  had  corresponding  MR  images  and 
mammograms  was  collected.  The  glandular  tissue  regions  in  the  MR  slices 
were  segmented  by  a  semi-automatic  method  and  the  percent  glandular  tissue 
volume  calculated.  Mammographic  breast  density  was  estimated  by  an 
automated  image  analysis  program.  It  was  found  that  the  correlation  between 
the  percent  dense  area  of  the  CC  and  MLO  views  and  the  percent  volumetric 
fibroglandular  tissue  on  MR  images  was  0.93  and  0.91,  respectively,  with  a 
mean  bias  of  4.4%.  The  high  correlation  indicates  the  usefulness  of 
mammographic  density  as  a  surrogate  for  breast  density  estimation. 


1.  Introduction 

Previous  studies  have  shown  that  there  is  a  strong  positive  correlation  between 
breast  parenchymal  density  on  mammograms  and  breast  cancer  risk  [1-3].  The 
relative  risk  is  estimated  to  be  about  4  to  6  times  higher  for  women  whose 
mammograms  have  parenchymal  densities  over  60%  of  the  breast  area,  as  compared 
to  women  with  less  than  5%  of  parenchymal  densities.  The  change  in  mammographic 
breast  density  is  therefore  often  used  as  an  indicator  for  monitoring  the  effects  of 
preventive  or  interventional  treatment  of  breast  cancer. 

Breast  cancer  risk  is  expected  to  depend  on  the  volume  of  glandular  tissue  in  the 
breast.  Mammographic  density  is  a  projection  of  the  volume  of  glandular  tissue  onto 
the  two-dimensional  image  plane.  To  better  understand  the  correlation  between 
mammographic  density  and  breast  cancer  risk,  it  is  important  to  investigate  the 
relationship  between  the  projected  areal  density  on  mammograms  and  the  volume  of 
glandular  tissue  in  the  breast.  In  this  study,  we  investigate  this  relationship  by 
analyzing  the  percent  volumetric  glandular  tissue  in  magnetic  resonance  (MR)  images 
and  the  percent  dense  area  in  corresponding  mammograms  for  the  same  breasts. 


2.  Materials  and  Methods 


Our  data  set  consisted  of  corresponding  MR  breast  images  and  mammograms  from 
37  patients  acquired  between  detection  and  biopsy.  The  MR  image  series  used  in  this 
study,  which  included  coronal  3D  Trweighted  pre-contrast  images  (coronal  sections 
2-5  mm  thick,  32  slices;  3D  Spoiled  Gradient-Recalled  Echo  (SPGR);  TE=3.3ms; 
TR=10ms,  Flip=40°,  matrix=256xl28,  FOV=28-32cm  right/left,  14-16cm 
superior/inferior,  scan  time=2  min  38  sec)  was  part  of  a  dynamic  breast  MRI  study. 
This  3D  SPGR  sequence  produced  full  volume  coverage  of  both  breasts  with 
contiguous  image  slices.  An  example  of  images  from  one  breast  is  shown  in  Fig.  1. 
Although  this  is  not  the  optimal  pulse  sequence  for  separating  water  and  fat,  the 
fibroglandular  parenchyma  (-water)  and  fatty  tissue  are  well  separated  with  this 
heavily  Trweighted  acquisition  and  therefore  the  series  was  chosen  for  this  study. 
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Fig.  1.  MR  images  of  the  right  breast  of  a  patient.  The  two-view  mammogram  of  the 
same  breast  is  shown  in  Fig.  2. 

We  have  developed  a  graphical  user  interface  that  displays  the  MR  series  and  the 
corresponding  mammogram  of  each  breast.  The  interface  allows  the  user  to  perform 
a  combination  of  manual  and  automatic  operations  to  segment  the  MR  images.  Each 
MR  slice  is  first  thresholded  to  separate  the  breast  from  the  surrounding  region.  For 
slices  close  to  the  chest  wall  where  no  clear  boundary  can  be  seen,  the  boundary  is 
manually  drawn  and  evaluated  by  radiologists.  Background  correction  [4]  using  the 
voxel  values  around  the  segmented  breast  region  is  employed  to  correct  for  the  non- 


uniformity  across  the  breast  area  due  to  the  breast  coil.  The  histogram  of  the  voxel 
values  in  the  breast  region  is  then  formed  and  interactive  thresholding  is  used  to 
segment  the  fibroglandular  tissue  from  the  fatty  tissue.  A  morphological  erosion 
operation  along  the  breast  boundary  then  excludes  the  skin  voxels  from  the 
calculation  of  the  fibroglandular  tissue  area  in  each  slice.  Finally,  an  integration  of 
the  fibroglandular  voxels  in  all  slices  relative  to  the  breast  volume  provides  the 
percent  volumetric  fibroglandular  tissue  in  the  breast. 

We  have  previously  developed  an  automated  image  analysis  tool  (Mammographic 
Density  ESTimator)  to  assist  radiologists  in  estimating  mammographic  breast  density 
[5].  MDEST  performs  dynamic  range  compression,  breast  boundary  tracking, 
pectoral  muscle  segmentation  for  the  MLO  view,  automatic  thresholding  based  on 
gray  level  histogram  analysis,  and  calculates  the  percent  dense  area  on  a 
mammogram.  We  found  that  the  correlation  between  the  computer-estimated  percent 
dense  area  and  radiologists’  manual  segmentation  was  0.94  and  0.91,  respectively,  for 
CC  and  MLO  views,  with  a  mean  bias  of  less  than  2%.  An  example  of  a 
mammogram  segmented  by  MDEST  is  shown  in  Fig.  2. 


Fig.  2.  Automated  breast  density  segmentation  from  mammograms.  Upper  row:  CC  view. 
Lower  row:  MLO  view. 


3.  Results 

Scatter  plots  of  the  percent  volumetric  fibroglandular  tissue  versus  the  percent 
dense  area  on  mammograms  are  shown  in  Fig.  3(a)  and  3(b)  for  the  CC-  and  MLO- 
view  mammograms,  respectively.  The  correlation  of  percent  dense  area  of  the  CC  and 
MLO  views  with  the  percent  volumetric  fibroglandular  tissue  on  MR  images  was 
found  to  be  0.93  and  0.91,  respectively,  with  a  mean  bias  of  4.4%. 


%  Area-CC  View  %  Area-MLO  View 


Fig.  3.  Correlation  of  %  volumetric  fibroglandular  tissue  on  MR  images  with  %  dense  area  on 
mammograms  for  37  patients.  The  left  and  right  breasts  are  plotted  as  separate  data 
points  on  each  graph.  The  dash  lines  are  linear  least  squares  fits  to  the  data  points. 


4.  Conclusion 

Mammographic  density  is  highly  correlated  with  the  volumetric  fibroglandular 
tissue  in  the  breast,  indicating  its  usefulness  as  a  surrogate  for  breast  density 
estimation.  The  computerized  image  analysis  tool,  MDEST,  is  useful  for  estimation  of 
mammographic  density.  The  automated  analysis  is  expected  to  contribute  to  the 
understanding  of  the  relationship  of  mammographic  density  to  breast  cancer  risk, 
detection,  and  prognosis,  and  to  the  prevention  and  treatment  of  breast  cancer. 
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